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Static packings of perfectly rigid particles with Coulomb friction are investigated theoretically and 
numerically. The problem of finding the contact forces in such packings is formulated mathemati- 
cally. Letting the values of the contact forces define a vector in a high-dimensional space enables us 
to consider the set of all possible contact forces as a region embedded in this same space. It is found 
that the boundary of the set is connected with the presence of sliding contacts, suggesting that a 
stable packing should not have more than 2M — 3iV sliding contacts in two dimensions, where M is 
the number of contacts and N is the number of particles. 

These results are used to analyze packings generated in different ways by either molecular dynam- 
ics or contact dynamics simulations. The dimension of the set of possible forces and the number 
of sliding contacts agrees with the theoretical expectations. The indeterminacy of each component 
of the contact forces are found, as well as an estimate for the diameter of the set of possible con- 
tact forces. We also show that contacts with high indeterminacy are located on force chains. The 
question of whether the simulation methods can represent a packing's memory of its formation is 
addressed. 

PACS numbers: 45.70.Cc 



I. INTRODUCTION 

The physics of granular materials involves two very dif- 
ferent length scales. The first length scale is associated 
with the size of the particles. If we wish to give an ac- 
curate value of the density, or describe the movement of 
the particles, it suffices to give the particle positions with 
an accuracy of some fraction of their radii. We call this 
length scale the "kinetic" length scale £kin, and take it 
to be of order the particle radius. When two particles 
touch, inter-particle forces at contacts are generated by 
tiny deformations that can be characterized by a second 
length scale, that we will call the "elastic" length scale 

One normally has l e \ <C 4ta' If one takes two marbles, 
£kin is about half a centimeter. But £ e i is not visible to 
the naked eye, as one can confirm by pushing the mar- 
bles together, and trying to observe their deformation. 
Because one often has £ e \ <C £kin, it is tempting to derive 
simplified numerical or theoretical approaches by taking 
the limit £ e i — > 0, corresponding to infinitely rigid par- 
ticles. One example is the inelastic hard sphere model, 
where collisions are assumed to be instantaneous. Instead 
of resolving the forces during a collision, one simply cal- 
culates the post-collisional velocities as a function of the 
pre-collisional ones. 

The inelastic hard sphere model opens the way to the 
application of kinetic theory and the use of event driven 
computer simulations. Both of these techniques have 
been applied successfully to a wide variety of granular 
flows pi 0, Q But the neglected length scale £ c \ takes 
its revenge in an unexpected way. If the collisions are 
dissipative, "inelastic collapse" can occur: there can be 
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an infinite number of collisions in finite time [fj. In 
event driven simulations, it is necessary to modify the 
simple inelastic hard sphere model to avoid this singular- 
ity. There are two general approaches. In the first ap- 
proach, each particle carries a clock that records the time 
of its last collision. When two particles collide, one checks 
their clocks to see if either one has had a collision less 
than some time t c ago. If so, the collision dissipates no 
energy, otherwise, the collision proceeds normally. The 
time t c corresponds to the duration of a collision [6j . Al- 
ternatively, one can make the restitution coefficient de- 
pend on impact velocity in such a way that the energy 
dissipation goes to zero as the impact velocity vanishes 
0. 

Another numerical method based on the approxima- 
tion Is! — > is contact dynamics (hereafter "CD") 
[H LLJJ ■ In this method, the contact forces are calcu- 
lated by requiring them to prevent particle interpenetra- 
tion and to minimize sliding. In this case, the neglected 
length scale takes its revenge by causing indeterminacy 
8J. In most cases, there are many possible force net- 
works that satisfy the constraints that are imposed on 
them. This raises several questions addressed in this pa- 
per: First of all, how big is the set of possible contact 
forces? Secondly, how are the forces chosen by CD dis- 
tinguished from all the other possible solutions? Finally, 
how do the forces chosen by CD differ from those cal- 
culated by soft-particle "molecular dynamics" (hereafter 
"MD") (ll). where the particle deformations are explic- 
itly treated? This paper addresses these questions. 

Another recently proposed a ppr oach similar to CD is 
the "force network ensemble" [T^l- The force network 
ensemble is the set of all possible force networks that 
could exist in a given configuration of rigid particles. 
In Ref. |l2j, this ensemble was sampled to obtain force 
distributions, which are compared with MD simulations. 
Parallels were drawn between the force network ensemble 
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and the ensembles of statistical mechanics, so that one 
could calculate properties of packings by averaging over 
the force network ensemble. But are all members of the 
ensemble equally likely to be realized? 

These questions have begun to be addressed. For ex- 
ample, it has been shown that the contact forces in static 
assemblies of frictionless grains can be uniquely deter- 
mined [lflj . In Ref. 0], the CD algorithm was adapted 
to sample the force network ensemble, allowing the au- 
thors to estimate its size. The authors found that the 
ensemble was not uniformly sampled, and that the force 
state generated by the dynamics had special properties. 
They also carried out a detailed study of the influence of 
tangential friction, and showed, that indeterminacy dis- 
appears in the limit of vanishing friction, consistent with 
Ref. [l3|. 

This paper takes a different, but complementary ap- 
proach. We investigate the structure of the force network 
ensemble mathematically, and show (in agreement with 
0|) that it is a convex set. In addition, we show that the 
boundaries of the set are associated with contacts where 
the Coulomb condition is marginally fulfilled. These find- 
ings place an upper limit on the number of such contacts 
that can exist in a static packing, and a lower limit on the 
dimension of the force network ensemble. We show how 
to locate the extremal points, enabling us to calculate 
the indeterminacy of the contact forces and to estimate 
the size of the force ensemble network. We also study the 
difference between the MD and CD calculation methods. 

This paper is organized into two main parts. Sec. [n] 
presents a mathematical formulation of the problem of 
finding the contact forces in a packing of infinitely rigid 
disks with Coulomb friction. It is shown that this prob- 
lem is equivalent to finding the intersection of a cone 
and a linear subspace in a high dimensional space. In 
Sec. IIIII we apply these ideas numerically to static pack- 
ings of about 100 particles, and answer several questions 
about the range of possible forces that could exist in the 
packing, and how they are related to the MD and CD 
solutions. 



II. MATHEMATICAL FORMULATION 
A. Definition of the contact matrix 

We study a system of N two-dimensional, circular 
grains at rest under gravity g in a rectangular container. 
We label each grain with a unique integer i, 1 < i < N. 
Particle i is characterized by its mass mi, radius r^, posi- 
tion fi, velocity Vi, momentum of inertia Ii, and angular 
velocity Ui. The fixed walls of the container could be 
considered as particles with infinite mass, but it is more 
convenient to simply leave them out of the analysis. 

Let M be the number of contacts between the N 
grains. Each contact can also be labeled with a unique 
integer a, 1 < a < M. Contact a is characterized by 
the two touching grains i and j. Given the positions of 



particles i and j, it is possible to define two unit vectors 
n a and t a that point in the directions normal to and tan- 
gent to the contact, respectively. At the contact, the two 
particles exert a normal force R a and a tangential force 
T a on each other. 

To calculate the motion of the particles, it is necessary 
to know the force fi and the torque T{ on each particle 
due to the contacts . Since fi and depend linearly on 
the contact forces, one can write 



cF. 



(1) 



Here, the contact forces and the forces experienced by the 
particles have been collected together into two column 
vectors F and f: 
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Note that f £ 



1>3N 



and F £ 



p2M 



The matrix c has 



dimensions 3iV x 2M, and is called the contact matrix. 
It contains information about which particles touch each 
other, and the geometry of the contacts. In Sec. Ill Bl we 
give explicit expressions for the components of c. Fol- 
lowing jSj, we call any particular value of F a "contact 
state" because it gives the state of all the contacts in the 
granular packing. 

Using this notation, we can easily write down the sys- 
tem of equations that must be solved in order to find the 
forces in a static granular packing under gravity. If the 
particles do not move, the contact forces must balance 
the gravitational acceleration: 



Mg. 



(3) 



Here, M £ R iN x W N is a diagonal matrix containing 
the masses and moments of inertia of the particles: 



/ mi 



M = 



\ 



mi 



h 



m N 



V 



m N 



(4) 



In I 



The vector g contains the gravitational accelerations of 
all the particles, organized in the same way as f in Eq. J2J). 

In contact dynamics, the unknowns are the contact 
forces F, not the forces on each particle, so we use Eq. 
to re-write Eq. ® as 



cF = Mg. 



(5) 
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If c were not singular, one could find a unique solution 
by inverting Eq. JSJ: 



F = -c _1 Mg. 



(6) 



But if c is singular, Eq. © has no unique solution, be- 
cause there exist vectors Fq ^ such that 



cF n 



0. 



(7) 



Physically, this means that there are contact states that 
exert no net force on the particles. This is the source 
of indeterminacy in granular packings j$\ . Once we have 
found a solution to Eq. J5j, we can construct an infinite 
number of solutions by adding multiples of Fo. How- 
ever, not all of these solutions are possible, for reasons 
discussed in Sec. Ill CI 

In general, c is singular, and its null space Co, plays 
a very important role in this paper. By considering the 
dimensions of c, one can establish a lower bound on the 
dimension of Co- c can be applied to any vector in R 2M , 
so its domain has dimension 2M. c maps this vector 
onto another vector in R 3Ar , so its range has dimension 
of at most 3iV. Since the dimension of the range and null 
space must add to the dimension of the domain, we have 



dimCn > 2M-37V. 



(8) 



If 2M < 3N, dim Co could vanish, but this corresponds 
to a coordination number of less than 3. Therefore, we 
expect that dim Co > 0. 



B. Construction of the contact matrix 

Suppose that particle i and j touch at contact a. Then 
one can define a unit vector n a normal to the particles' 
surfaces at the contact: 



(9) 



To discuss tangential forces, we must define a unit tan- 
gential vector such that h a ■ t a = 0. Given h a , there are 
two possibilities for t a , but t a can be uniquely defined by 
imagining that the two dimensional space is embedded in 
three dimensions, and writing 



(10) 



where z is the unit vector, pointing upwards, perpendic- 
ular to the two-dimensional plane. 

With these definitions, it is now possible to write the 
forces A/i Q and Afj a due to contact a on particles i and 
j- 

^fia — ^a^a T(yt ai Afj a R a ii a T a t ai (11) 

But this notation is awkward, because it is necessary 
to distinguish the two particles of the contact. One of 



the particles is "first" (particle i) and the other is "sec- 
ond" (particle j). The choice of which particle is first 
is arbitrary, but once the choice is made, it must not 
be changed. Accordingly, we introduce the symbol Xia 
defined by 



X 



1 if particle i is first in contact a, 
— 1 if particle i is second in contact a, 
if particle i does not participate in contact a. 

(12) 

For each contact between two grains, one element of \ 
is 1, and another is —1. Contacts between a wall and a 
grain contribute only one nonzero element to x- Using 
the x symbol, Eqs. Ijllfl as a single equation: 



A/fca = Xka{R a n a + ' 



(13) 



This equation holds for 1 < k < M, so the total force on 
a particle can be written as a sum over all the contacts: 



fk — 2, Xka(R a n a + T a t a ) 



a=l 



(14) 



This equation can be cast in the form of a matrix multi- 
plication. From Eq. (|14|) and an analogous equation for 
the torques, it is possible to deduce the components of c. 
c is a N x M matrix of submatrices Ci a , where 

XioJ^ax Xict^a 

Xicn^ay Xiatoty \ ' (-^^) 
\Xia\n 



C. The contact conditions 

Eq. © does not give a complete description of motion- 
less granular packings. Granular packings are nonlinear 
because only certain contact forces are physically possi- 
ble. For dry granular materials with Coulomb friction, 
two conditions need to be met: 



R a > 0, and \T a \ < [iR a , 



(16) 



for a — 1,...M. The first condition says that there 
are no attractive forces, only repulsive ones. The second 
condition states that the tangential force cannot exceed 
fi times the normal force, where the constant fj, is the 
Coulomb friction ratio. Let us define K to be the set 
of all contact states satisfying Eq. (jlf)|l . In Fig. ^ we 
show the cross section of K, cut by the (R, T) plane of a 
contact. 

Later in this paper, our calculations will be greatly 
simplified because K is a convex set. This can be seen 
from Fig.^ since K is a cone with its vertex at the origin. 
The convexity of K can also established by a simple proof. 
Consider two points Fa, Fb el. The points 



XF A + (1 - X)F B , 0<A<1. 



(17) 
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FIG. 1: Cross section of K cut by the (R,T) of a contact. 
Contacts that satisfy Eq. I|16|l must lie in the shaded triangu- 
lar region. 

lie on a straight line between and Fb- One can show 
that Fa £ K as well, showing that K is convex. 

As we shall see, contacts where an equality holds in 
Eq. JTBJl, i.e., where R a — or \T a \ — fiR a , play a spe- 
cial role in limiting the indeterminacy. Contacts with 
R a = are called "non-transmitting" contacts because 
the particles touch, but exert no force on each other. 
Contacts where \T a \ = fiR a holds will be called "sliding 
contacts" even if there is no relative motion. This term is 
used because this equality holds when two particles slide 
relative to one another. 



D. The set of possible contact states 

We are now in a position to find the contact forces in 
a static granular packing and understand how indeter- 
minacy arises. Given an arrangement of particles and 
external forces acting on each particle (e.g. gravity), one 
can calculate the contact matrix c and the vector Mg. 
Then one first searches for a "particular solution" Fi such 
that cFi = Mg. These forces are necessary to cancel 
the external forces. The particular solution Fi can be 
made unique by requiring that it be orthogonal to every 
vector in Co- 

In general, Fi will not obey the contact conditions 
Eq. I|lrj|) . so one must find some Fo G Co such that 

Fo + FieK. (18) 

The sum F + Fi is a possible contact state. There could 
be many vectors F G Co that satisfy Eq. (|T%|) . so the 
solution may not be unique. On the other hand, not all 
vectors from Co will satisfy Eq. (|18|) . 

Let F be the set of all such contact states F = Fo + Fi 
satisfying Eq. (|18|l . Then F is the set of all contact states 
that could be observed in a given granular packing. We 
can summarize the requirements of a contact state by 
writing 

F = |f + Fi | F + Fi 6 K, F € C , cFj = Mg, 



Fi • f£° = for i = 1. . .dimCo j.(19) 

Here jF^\i = 1 . . .dim Co} denotes a basis of Co- In 
Ref . [12| , F is called the "force network ensemble" . 

It is clear that the set F is also convex. Let F a and F& 
be members of F. Then the intermediate points Fa have 
the form 

Fa = AF a + (1 - A)F 6 = Fr + AF , a + (1 - A)F 0)6 
= F!+F ,a, (20) 

where we have used the decompositions F a = Fi + Fo, a 
and Fb = Fi + Fo.&. It can be shown that Fa satisfies all 
the conditions in Eq. (|19(l : Fa S K because K is convex 
and Fo.a £ Co because Co is a linear subspace. The 
convexity of F was first noted in Ref. |15| , and shown in 
Ref. 

If F is empty, the granular packing is unstable, and the 
particles will move. If F contains only one point, there 
is a unique contact state, and there is no indeterminacy. 
Finally, F can contain many points. In this case, there is 
no unique solution. 

E. The Boundary of F 

The boundary of F is analogous to the "yield surface" 
in the elastoplasticity theory of soils ^tJ- K a system 
crosses the boundary, and leaves F, the packing is no 
longer stable, and the particles will move. When this 
happens, a new contact matrix must be constructed, in- 
cluding perhaps new contacts, and F and Co will change. 
For our purposes, we are interested in the boundary of F 
because it contains the extremal points, where the con- 
tact forces are maximized or minimized. 

To investigate the structure of F, we will use the con- 
cept of an n dimensional neighborhood of point. We say 
a point FgF has an n dimensional neighborhood if there 
exist at most n linearly independent vectors V such that 
one can find a mlu < and a max > satisfying 

F + aV e F, for all a min < a < a max . (21) 

As an example, consider a line segment embedded in two 
dimensional space. Each point on the line segment has a 
1 dimensional neighborhood, and the end points have 
dimensional neighborhoods. 

Now let us apply this concept to our set F. Let us 
suppose that that F contains a contact state Po with 
no sliding or non-transmitting contacts. A multiple of 
any vector in Co can be added to Po, as long as it is 
small enough. Therefore, there are n — dim Co linearly 
independent vectors V satisfying Eq. 1|21[1 . meaning that 
Po has an n dimensional neighborhood. 

Next suppose that we have a point Pi £ F with exactly 
one sliding contact. Let us label that sliding contact 
and assume that we have fiR^ = Points P a = 

Pi + aV obey 

M< - T« = MflW - T« + a {id" ~ if >) , (22) 
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where R^\ R@~\ and R^ are the appropriate compo- 
nents of P a , Pi, and V respectively. Since contact f3 is 
sliding in the state Pi, we have fJ,Rg — — 0. Now, 
note that a takes on both positive and negative values. 
If P a is to satisfy the contact conditions, we also need 
— Ti, = 0, i.e., contact (3 must be sliding in V 
also. Thus the sliding contact puts a constraint on the 
vectors V that can be used. However, given any two 
vectors from Co, one can construct a linear combination 
satisfying this constraint. In this way, n — 1 linearly in- 
dependent vectors can be constructed, so Pi has an n — 1 
dimensional neighborhood. Similar reasoning can be ex- 
tended to show that a contact state with two sliding con- 
tacts has ann-2 dimensional neighborhood, and finally, 
a contact state with n sliding contacts is a dimensional 
neighborhood, that is, it is an extremal point. 

So far, non-transmitting contacts have not been con- 
sidered. But it is easy to incorporate them, because 
they can be considered as two sliding contacts superim- 
posed on each other, i.e., the contact obeys T = fiR and 
T = —fiR at the same time. 



F. Quantifying indeterminacy 

The indeterminacy of the granular packing is deter- 
mined by the size and shape of the set F. Since F is 
a convex set with a finite number of extremal points, 
one possible approach would be to locate all its extremal 
points, but the large number of such points makes this 
unfeasible. Therefore, we adopt an alternative approach. 
We locate the subset F ext of extremal points where a 
component of some contact force attains its maximum or 
minimum possible value. Specifically, for each contact a, 
four different extremal points are found: the state where 
R a is maximum, then where R a is minimized, and then 
the states where T a is maximized and then minimized. 
(Note that minimizing a normal force means making it 
approach as closely as possible, but minimizing a tan- 
gential force means making it approach -co.) 

Once F ext has obtained, different measures of indeter- 
minacy can be extracted. One possibility is to calculated 
the range of possible forces that a given contact a can 
take on: 

c ^max.Q -^min,a c -Mnax,a T mma 

OR.a = ZZ ! OT,a = ZZ , 

mg mg 

(23) 

where i? ma x,a and i? m in iQ are the maximum and mini- 
mum possible values of the normal contact force at con- 
tact a. To obtain a dimensionless number, we divide by 
mg, the average weight of a particle. We call Sr and 
St the "local indeterminacies" because they quantify the 
ambiguity of the force at one contact. 

One contact force cannot be maximized independently 
of the others. One could therefore ask how much the 
entire contact force state must change when we bring one 
contact force from its minimum value to its maximum. 



As an alternative to Sr and St, one could calculate 

j 1 1 Pi?,, max, a Fi?,min,a|| 



j ||FT,max,a PT,min,a 1 1 (c\a\ 
«T,a = ZZ ■ (24 ) 

mg 

Here, F_R.max,a is the contact state where the normal 
force is maximized at contact a, and the other contact 
states in Eq. I|24j) are defined analogously. These quan- 
tities estimate the "diameter" of F. Note that F is em- 
bedded in a 2M dimensional space, so that 2M different 
"diameters" can be calculated. We call da and d,T the 
"global indeterminacies" . 

G. Algorithm 

We now give the algorithm used to find the maximum 
and minimum possible forces at a given contact. The 
first step is to locate an extremal point. Let us suppose 
that we begin with a point with an n-dimensional neigh- 
borhood in F. Given this point, an extremal point can 
be found in the following way: Pick one of the n vectors 
from the basis of Co- Then, starting from the interior 
point, move in that direction until a sliding contact is 
detected. Then n — 1 linearly independent vectors can be 
constructed out of the basis of Co, all of which preserve 
the status of the sliding contact. Pick one of these vec- 
tors, and proceed in its direction until a second sliding 
contact is detected (or the first sliding contact becomes 
non-transmitting) . Then n — 2 linearly independent vec- 
tors can be built which obey these two constraints. Con- 
tinuing in this way, we will eventually reach a point where 
there are n constraints, and no vectors can be constructed 
that respect all of them. This is an extremal point. 

Once the extremal point has been reached, its neigh- 
boring extremal points can each be identified. Recall 
that an extremal point is characterized by n constraints 
arising from n sliding contacts. If we relax one of these 
constraints, there is one direction in Co that respects all 
the other n — 1 constraints. If we move away from our 
extremal point in this direction, we will eventually en- 
counter a new extremal point. This is a neighboring ex- 
tremal point, connected to the current point by an edge. 
Since there are n possible constraints to relax, each ex- 
tremal point will have n neighbors. Each of these neigh- 
bors can be checked. If none of them are "better" than 
the current point (in the sense that the relevant contact 
force is greater or less), then the current point is the best 
point. If any one of the neighboring points is better, we 
move there and repeat the process. The convex struc- 
ture of F guarantees that there are no local minima or 
maxima that would trap the algorithm. 

The algorithm must deal with a number of practical 
difficulties. For example, it can happen that a contact 
must always be sliding or non-transmitting. It is neces- 
sary to detect this, because such a situation reduces the 
dimension of F. Therefore before beginning to search for 
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FIG. 2: The particle i is supported by two walls through 
contacts a and j3. Gravity pulls the particle downwards. The 
angle <j> suffices to characterize the geometry of this simple 
granular packing. 



extremal points, we first try to locate a point without 
sliding contacts. This can be difficult, because the simu- 
lations often yield points with many sliding contacts. But 
given such a point, one can construct a vector belonging 
to Co that preserves all the sliding contacts except one, 
and then moving along this vector so that the number 
of sliding contacts is reduced by one. This process can 
be repeated. Sometimes this does not work, because the 
simulation yields a point in a tight, multi-dimensional 
corner. In this case, the contact dynamics iterative solver 
can be used to generate an alternative starting point. 

In about 1% of the cases analyzed in Sec lIIII ex- 
tremal states are found with a different number of slid- 
ing contacts than expected. This may happen when one 
constructs linear combinations that satisfy a given con- 
straint, and by chance, satisfy several other constraints 
at the same time. Another possibility is that it is diffi- 
cult to maintain sufficient numerical accuracy during the 
construction of the linear combinations. Recall that one 
cannot test strict equalities with floating-point numbers; 
one should always check that equality conditions such as 
T = ±/ii? are satisfied to a certain tolerance. This may 
cause an occasional overestimate of the number of slid- 
ing contacts. However, since this situation occurs only 
rarely, it does not affect the conclusions of this paper. 



H. Two Contacts 

As a simple application of these ideas, consider the 
granular packing in Fig[2 where M — 2 and N = 1. The 
equations of static equilibrium are 



R a sin 4> + T a cos 4> — Rp sin 4> + Tp cos <fi 
R a cos 4> — T a sin <j> + Rp cos <j> + Tp sin <f> 

rT a + rTp 

Comparing this to Eq. (J3J, we have 




0, 

mg, 

0. (25) 



(26) 



This matrix has a null space that has at least one dimen- 
sion. Let us find Fo by solving the system of equations 
cF = 0. The result is 



'o = 



V 




(27) 



Note that Fo corresponds to the horizontal components 
of the contact forces canceling each other. 

The particular solution can be found by solving cFi = 
Mg, and then requiring that Fq • Fi = 0. The result is 



/ cos< 



mg 
2 



— sm< 

COS0 

smd> 



(28) 



which simply expresses the requirement that the verti- 
cal component of the contact forces cancel gravity. The 
contact states in F all have the form 



Fi +aF . 



(29) 



Applying the contact conditions puts restrictions on the 
values of a which are allowed. For example, requiring 
Ra > and Rp > means that a must satisfy the in- 
equality 



mg 

a > cot (p. 

And the Coulomb condition becomes 

mg . 



mg 
\i— cost 



/ia sin i 



< 



a cos < 



(30) 



(31) 



Working out the various cases connected with the abso- 
lute values, one obtains a lower bound for a. 



a > a n 



mg 
2 



fi — tan (f, 
1 + /i tan < 



(32) 



If this condition is fulfilled, than Eq. I|3(J[) is always satis- 
fied as well. When tan0 < l/fi, there is an upper bound 
for a: 



a < a E 



mg 



/j + tan cj) 
1 — fi tan < 



(33) 



If this condition is fulfilled, then Eq. H3Ufl is satisfied as 
well. When tan0 > there is no upper bound on a; 
a can be arbitrarily large. 

Now the indeterminacy of this packing can be cal- 
culated. Selecting the appropriate components from 
Eq. we have 



R a — R 



■0 



mg 



(cos (j) 



(34) 



The maximum and minimum possible forces can be ob- 
tained by setting a equal to a max or a m ; n respectively. 
When this is done, one obtains 



Rn 



mg 



mg 
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Grains generated 
on a grid. 
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FIG. 3: The measures of indeterminacy Sr (solid line) and St 
(dotted line) . 4> is the angle shown in Fig. [5] and n was taken 
to be 0.3. At (j> = 0, <5i = fi and 82 = Both measures 
diverge at tan</> = or « 0.417T. 



CD-2 



MD-2 



~~ 2 J. 9 * 2 i' v oiJ / 
cos z <p — sin 

In the same way, one has 

— T Q = Tjg = (sin <f> ~ a cos 0) . (36) 
Inserting a — a max and a — a m ; n , we obtain 

— (37) 
cos z — //^ sin <p 

The behavior of Ecis. H35|) and l|37[) are shown in Fig. EJ 

III. NUMERICAL APPLICATION 

A. Overview 

We have investigated numerically the indeterminacy of 
granular packings with N = 95 particles. Fig. 0] shows 
the procedure used to generate the various configurations 
considered in this section. First, N = 95 grains were 
placed on a grid inside a rectangular box of size L x x L y . 
To prevent the formation of regular arrays of particles, 
the radii of the particles are uniformly distributed in the 
interval [0.7r max , r m ax]- Then, the grains were allowed to 
fall, under the influence of gravity, and form a packing 
at the bottom of the box. This process was simulated by 
both CD and MD, yielding the CD-I and MD-1 config- 
urations. Both simulations continue until the kinetic en- 
ergy decreases to a negligible value, or the elapsed time 
reaches 6y/L y /g, i.e. about twice the time a particle 
needs to fall a distance L y . (This second condition is 
needed because occasionally a particle falls to the bot- 
tom and starts to roll with a small amount of kinetic 
energy that is enough to violate the first condition, but 



FIG. 4: A diagram showing the relationship between the var- 
ious configurations analyzed in this section. First, 95 grains 
were placed on a rectangular grid. Then, the grains were al- 
lowed to fall, under gravity, and settle at the bottom of the 
container. This process was simulated by both CD and MD, 
yielding the CD-I and MD-1 configurations. Then the CD-I 
configuration was taken as the initial condition for two more 
simulations, yielding CD-2 and MD-2. 



not large enough so that it reaches another particle or 
a wall in a reasonable amount of time.) Then the CD-I 
configuration is used as the starting point for two more 
simulations, yielding the CD-2 and MD-2 configuration, 
according to the simulation method used. In this second 
CD simulation, the initial guess for the contact forces is 
F = 0. This procedure was repeated sixty times, yield- 
ing a set of 240 configurations. The 60 initial conditions 
are distinguished by choosing different particle radii each 
time. 

In the MD simulations, particles interact via linear, 
damped springs in both the normal and tangential direc- 
tions. When a contact becomes sliding, it is assumed that 
the tangential spring remains stretched at its maximum 
length. The springs have stiffness 1.2 x 10 5 rng/r and 
damping constant 7.7 x 10 4 (Wig / V) y/rjg ■ This choice of 
parameters leads to an average particle overlap of 7x 10 -6 
of a particle radius. Usually particles are much softer in 
MD simulations, but extremely hard particle were used 
here to approach as closely as possible the CD simula- 
tions (average overlap: 6 x 10~ 8 of a particle radius). 
Due to the hardness of the particles, the MD simulations 
were very slow, taking roughly 100 times as long as the 
CD ones. In the CD simulations, the normal and tan- 
gential resitution coefficients were set to 0. In all cases, 
the Coulomb friction ratio was /j, = 0.3. 

Then the resulting configurations are analyzed. The 
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FIG. 5: Three different contact states for a configuration with 
N = 95 particles. The CD-I state obtained from contact dy- 
namics is on the left. There are M = 160 contacts. Then 
this configuration was allowed to relax in a molecular dynam- 
ics simulation, resulting in the MD-2 configuration shown in 
the middle. During the relaxation, the number of contacts 
increases to M = 166. The state on the right is the extremal 
state of the CD-I configuration with the maximum norm. The 
thickness of the lines connecting the centers is proportional 
to the normal force, and tangential forces are shown by lines 
tangent to the particle surfaces. 



matrix c was constructed, and a singular value decom- 
position was used to extract its null space. Then the set 
F cxt defined in Scc. lIIF| is found using the algorithm pre- 
sented in Sec. Ill Gl Some examples of the contact states 
found are shown in Fig. EI The left hand panel is the CD- 
1 state, the middle panel is the MD-2 state, right hand 
panel shows one of the elements of F cxt - One can see 
already that all three states are different, but the CD-I 
and MD-2 states are closer together than either one with 
the extremal state. 

Fig. shows the dimension of F as a function of the 
number of contacts. All points fall onto or just above 
the line dimF = 2M — 3N, confirming the prediction 
that dimF > 2M - 3N. Fig. also shows that dimF 
never exceeds 2M — 3N by more than 3. This means 
that it is reasonable to use dim F w 2M — 3iV to estimate 
dimF. There is also no difference between the classes of 
configurations, except that MD simulations tend to have 
more contacts than the CD ones. 

Fig. m> shows the number of sliding contacts observed 
in the different configurations. A contact is considered 
sliding if fiR - \T\ < emg. For Fig. Eb, e = 10 -9 . Non- 
transmitting contacts (R, \T\ < emg) are counted as two 
sliding contacts. The number of sliding contacts is al- 
ways less than or equal to 2M — 37V, consistent with our 
analysis of the boundary of F in Sec. Ill El The difference 
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FIG. 6: (a) the dimension of the set of possible contact states 
F as a function of the number of contacts M. The solid line 
shows the lower limit dimF = 2M—3N. (b) the number M a of 
sliding contacts as a function of dimF. The straight line shows 
M s = 2M - 3N. A contact is considered sliding if (j,R- \T\ < 
(-Trig, with e = 10 -9 . Non-transmitting contacts (R, T < emg) 
count as two sliding contacts, consistent with our discussion 
in Sec. Ill El Four different families of simulations are shown: 
triangles - MD-1, circles - MD-2, crosses - CD-I, x's - CD-2. 



between the different classes of configurations becomes 
clear. The CD simulations have very few sliding con- 
tacts. The two classes of MD simulations are also well 
separated from each other, with the MD-2 configurations 
having the most sliding contacts. This difference between 
the two configurations shows that the MD simulation is 
able to retain a memory of how it was formed. The MD- 
1 simulation was generated by letting the particles fall 
from a given height, whereas the MD-2 simulation was 
started with the particles almost in their final position. 
Therefore, much more energy was dissipated during the 
MD-1 simulation than the MD-2 simulations. The differ- 
ence in the number of sliding contacts is a sign of their 
different histories. 

Note that the CD simulations lack this kind of mem- 
ory. A closer examination of the CD-2 simulations re- 
veals that there are initially many sliding contacts, but 
this number decreases rapidly to values nearly zero, as 
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FIG. 7: Histograms of local indeterminacy Sr (thick lines) and 
St (thin Unes), denned in Eq. l]2.'-j> . for the CD configurations 
(solid lines), MD configurations (dashed lines), plotted semi- 
logarithmically. The thick lines show 5r and the thin lines 
St- The distributions are normalized so that their integral is 
unity, hence we call them probability density functions. 

shown in Fig. [SJs. There is one exceptional simulation, 
indicated by the cross near M s « 20, dimF « 30, where 
the system seems to be trapped in some corner of F and 
unable to escape. (The nearby CD-I simulation prob- 
ably represents a similar situation, but the two points 
were generated by different random number seeds, so it 
is probably a coincidence that they are so close together.) 



B. Measurement of Indeterminacy 

1. Local indeterminacy 

The local indeterminacy presented in Sec. Ill Fl can be 
calculated. In Fig.[7J we show the distribution of Sr and 
St for the four different families of configurations. In all 
cases, the distributions are exponential, except where a 
sharp peak appears near 5r, 8? = 0. Note that similar 
exponential tails are observed in contact force distribu- 
tions. 

The MD distributions show slightly larger indetermi- 
nacies than the CD ones, independent of the method of 
generating the configuration. No contacts were observed 
with infinite indeterminacies, although it is theoretically 
possible, as was shown in Sec. Ill HI Infinite indetermi- 
nacy may exist in only special packings with very small 
number of particles. 



2. Global indeterminacy 

Let us now consider global indeterminacy. Histograms 
of the d,R and dx, defined in Eq. Q24[l. are shown in Fig. El 
for the three different families of configurations. This 
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FIG. 8: Histograms of global indeterminacy da (thick lines) 
and dr (thin lines), defined in Eq. I|24|l . for the CD simulations 
(solid lines), and the MD simulations (dashed lines). 

measure of indeterminacy has very different properties 
from the previous one. In Fig. [7J the most probable val- 
ues of the indeterminacy were small, but in Fig. [H] the 
probability density function presents two maxima, one 
close to dn,dT — 0, and the another well separated from 
the smallest values. The second maximum is by far the 
largest, so maximizing or minimizing a contact force usu- 
ally involves changing many forces throughout the pack- 
ing, even when the change at the contact in question 
is small. The maxima in Fig. [5] can be taken as crude 
estimates of the diameter of F, indicating that F has a 
diameter of approximately 60mg in the CD simulations 
and 75mg in the MD simulations. Consistent with Fig.JTJ 
the MD simulations show slightly higher indeterminacy 
than the CD ones. Finally, note that the curves for the 
normal and tangential components are nearly the same 
in Fig. [SI but clearly different in Fig. [7J 



3. Force Chains 

One of the most remarkable characteristics of granular 
packings is that most of the force is carried by a small 
fraction of the contacts, which are organized in linear 
structures called force chains. Examples of force chains 
are clearly visible in Fig. |SJ This figure also suggests that 
indeterminacy is also concentrated along force chains - 
it is the force bearing contacts which differ the most be- 
tween the three panels. We can confirm this impression 
by dividing the contacts into two different classes, those 
with above average forces and those with below average 
forces. In Ref. [lq. it was shown that these two popu- 
lations of contacts have different properties, the former 
being associated with force chains. (In our case, these 
two classes do not have exactly the same meaning as in 
Ref. 0, as the stress is not uniform throughout the 
packing. Nevertheless, it is still possible for us to iso- 
late the contacts in the force chains, at least in the lower 
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FIG. 9: Probability distributions of the local indeterminacies 
Sr (solid lines) and St (dashed lines) for contacts with R < R 
(thin lines), and for R > R (thick lines), where R = 9.3m<? is 
the average normal force. Only data from the CD-I simula- 
tions are shown; the others yield similar curves. 
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4- Alternative measurements of indeterminacy 

In Ref. indeterminacy is measured in a different 
way. A configuration is given to the CD simulation pro- 
gram, and iterative solver of the CD simulation is asked 
for a possible solution to the forces. This solution is then 
perturbed slightly, and given to the iterative solver as 
an initial guess, and a new solution is obtained. This is 
repeated many times, and the force network ensemble is 
sampled in the same way as different statistical mechan- 
ical ensembles are sampled in Monte Carlo simulations. 
Of course, there is no guarantee that that this method 
will weight appropriately the different regions of F or even 
that it will explore all parts of F. 

We have carried out this procedure on our packings. 
For each configuration, we have obtained 500 different 
solutions. The solutions are perturbed by multiplying all 
the contact forces by a random number uniformly dis- 
tributed between 0.5 and 1.5. Then the center F of F 
can be estimated by averaging over all 500 points, and 
the radius r of F can be estimated by the variance: 

/ 500 \ 1 / 2 

-(^X>-F°] 2 J ■ (38) 

Here, Fj are the solutions obtained from the iterative 
solver. Averaging over all 60 configurations in each class, 
we obtain r = (33 ± 1.4)mg for the CD-I configurations, 
r = (34 ± l.3)mg for CD-2, r = (41 ± l.7)mg for MD-1, 
and r = (40 ± 1.5)mg for MD-2. These values are quite 
close to half the diameter of F estimated from Fig. [8] 
In all cases, the standard deviation of r is about 12mg, 
which is also consistent with Fig. |H1 This suggests that 
one can indeed sample F in this way. 



FIG. 10: Probability distributions of the global indetermina- 
cies da (solid lines) and dr (dashed lines) for contacts with 
R < R (thin line), and for R> R (thick line), where R is the 
average normal force. Only data from the CD-I simulations 
are shown; the others yield similar curves. 



part of the packing.) In Fig. we show the distribu- 
tions of local indeterminacy in the CD configurations for 
each class of contact. The two classes yield quite differ- 
ent distributions. The peak at Sr,St = is due entirely 
to the contacts with below average force, while the con- 
tacts with large indeterminacy have above average forces. 
Thus, force chains are also "indeterminacy chains" . 

We can examine the distribution of global indetermi- 
nacy as well. This is done in Fig. 1 101 The small maximum 
near d,R,dT = is due entirely to contacts with R < R. 
The main maximum has contributions from both classes 
of contacts, but removing the weak contacts causes the 
maximum to shift towards larger values. 



C. Comparison between MD and CD 

1. Distance between the states 

Next, we would like to examine more closely the differ- 
ence between the forces calculated by MD and CD. One 
way to do this is to compare the CD-I, CD-2, and MD-2 
configurations, where the particles positions are nearly 
identical. The distance between two states Fa and Fb 
is simply \\Fa — Fb||. If Fa and F# have different num- 
bers of contacts, one can insert 0's into the appropriate 
places so that they have the same dimension. 

We find that the distance between CD-I and CD-2 is 
(50±3)m5, while (75±A)mg separates CD-I and MD-2. 
Thus, the MD moves the forces farther away from the ini- 
tial state than CD does. But the MD-2 and CD-2 states 
are separated by only (45 ± 3)mg, indicating that they 
move in approximately the same direction. Note that all 
these distances are less than, but of approximately the 
same magnitude as, the diameter of F. Hence changing 
from CD to MD or erasing the memory of CD causes 
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R=(R-R )/(R -R . ) 
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T =(T -T . )/(T -T . ) 
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FIG. 11: Histograms of R* and T*, defined in Eq. H39H . for 
the four different configurations. The solid lines show CD 
simulations, and the dotted lines show MD simulations. The 
thick lines show the results of dropping the particles (CD-I 
and MD-1), and the thin lines show the result of obtaining the 
forces with very little particle movement (CD-2 and MD-2). 



a perturbation in the forces of roughly the same size as 
their indeterminacy. 
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FIG. 12: Histogram of the normal force 7? in the simulated 
states and various classes of extremal states. Sixty configu- 
rations of N — 95 particles were considered. The curves for 
the extremal states have less noise because each configuration 
contributes M extremal states, but only one contact dynamics 
state. 



of configurations. Surprisingly, these distributions de- 
pend more on the history of the configuration than on 
the simulation method. When the particles are dropped 
and allowed to settle (CD-I and MD-1), the normal forces 
are larger, relative to their minimum and maximum pos- 
sible values, than when the particles are simply placed 
into their final positions (CD-2 and MD-2). In the tan- 
gential case, the CD-I and MD-1 states show more val- 
ues clustered around the middle of the allowed interval 
than the CD-2 and MD-2 states. On the other hand, the 
CD-2 and MD-2 states have more contacts at their max- 
imum or minimum values. These results show that the 
CD method is capable of representing the history of the 
packing, even though this is not reflected in the number 
of sliding contacts. 



D. Contact Force distributions 



2. Relation to extremal states 

To see more precisely where the state found by the 
simulation stands in relation to the extremal states, let 
us consider the following quantities: 



R* 



Rs 



Rn 



R„ 



T., 



(39) 



where R s im and T s - lm are the contact forces found by the 
simulation, and i? ma x and T max are the maximum values 
these forces could attain in this configuration, while i? m in 
and T m j n are the minimum forces. i?» = means that 
the simulation chose the minimum possible force while 
R* = 1 means that it chose the maximum. In Fig. 1111 
we show histograms of i?* and T* for the four families 



1. Extremal vs Simulated states 

The contact force distributions in the various types 
of contacts states are compared in Fig. ^] The ex- 
tremal and contact dynamics distributions coincide for 
R < SOrrig, but then separate, with the extremal states 
exhibiting more contacts with large forces. This is true, 
even of extremal states found while minimizing R. All 
the distributions are approximately exponential, similar 
to those found in other studies [13, Ha] • This result sug- 
gests that the exponential contact force distributions are 
a property of all members of F. This means that these 
exponential tails are probably due to some property of 
Eq. |JS}, and can be studied using the force network ap- 
proach. 
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FIG. 13: Histogram of the normal force R for the different 
families of simulations: CD-I - thick solid line, CD-2 - thin 
solid line, MD-1 - thick dashed line, MD-2 - thin dashed line. 

2. MD vs CD 

Finally, we compare in Fig. ll3l the normal contact force 
distributions for the four different families of simulations. 
Both CD-I and MD-1 yield similar curves, obeying an 
exponential fall off out to the largest observed values of 
R. On the other hand, CD-2 and MD-2 have fewer con- 
tacts at these large values. This suggests that during the 
formation of the packing, some of the kinetic energy is 
stored as elastic energy in the force chains. If the packing 
is formed with very little kinetic energy (as in the case of 
CD-2 and MD-2), fewer large contact forces are present. 

IV. CONCLUSIONS 

One conclusion that can be drawn from this work is the 
importance of sliding contacts. In Sec. Ill Fl we showed 
that they are associated with the boundary of F, and 
hence are a sign that the packing is close to yielding. 
Furthermore, our results suggest that there should be 
fewer than 2M — 3N sliding contacts (counting non- 
transmitting contacts as two sliding contacts) in a sta- 
ble packing, and no counterexamples were found among 
the 240 configurations that were examined. If a granular 
packing is slowly loaded, our work predicts that the num- 



ber of sliding contacts will increase, and reach 2M — 37V 
when the packing yields. The MD algorithm produces 
many more sliding contacts than CD. This suggests that 
packings under CD are much more stable than under MD. 

We were able to calculate the local indeterminacy, that 
is the range of values a given contact force can assume. 
We found that the contacts with large indeterminacy are 
also those contacts that make up force chains. There- 
fore, the primary origin of indeterminacy is that the am- 
plitudes of the force chains can change. This result also 
suggests that force chains could be understood by inves- 
tigating the null space of the contact matrix c, since the 
difference between any two allowed states belongs to this 
set. 

The global indeterminacy measures how much the en- 
tire network must be adjusted in order to maximize or 
minimize a force at a given contact. The reorganizations 
required for most contacts are significant, even when the 
local indeterminacy is small. The global indeterminacy 
can also be used to estimate the diameter of F. This di- 
ameter in turn allows one to appreciate the magnitude of 
changes occurring within the force network. For exam- 
ple, we saw that changing the simulation method changes 
the force network by an amount roughly equivalent to 
the diameter of F. Erasing the memory of a CD simu- 
lation changes the forces by roughly half as much. Our 
estimates of the diameter of F are consistent with those 
obtained using a Monte Carlo-like procedure to sample 
F. 

Finally, we made several observations about how a 
packing's "memory" is formed. When a packing is formed 
violently, with much kinetic energy, some of this energy 
ends up stored in the contacts. Such packings exhibit 
stronger force chains and larger contact forces than pack- 
ings formed gently, with very little kinetic energy Both 
CD and MD simulations show this effect, although it is 
more significant in the MD simulations. 
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